1 Fitting Archetypes to ROSMAP RNAseq data
Annat Haber
2020-09-09
This script conducts archetpye analysis following
https://vitkl.github.io/ParetoTI/
https://royalsocietypublishing.org/doi/10.1098/rstb.2017.0105
https://www.nature.com/articles/nmeth.3254
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
knitr::opts_chunk$set(cache=TRUE, cache.extra = R.version, warning=FALSE, message=FALSE)
reticulate::use_condaenv("reticulate_PCHA", conda = "auto", required = TRUE)suppressPackageStartupMessages({
library(ParetoTI)
library(ggfortify)
library(GGally)
library(cowplot)
library(Matrix)
library(tidyverse)
})Preparing data and metadata
Loading data downloaded directly from synapse, as well as Nikhil’s subtypes, and merging it all together.
Diagnosis variable combines cogdx, braak, and cedar scores following syn8456629
Dementia that isn’t MCI or AD (cogdx=6) is removed so that “Other” includes only MCI.
Data is arranged by rnaseq_id (rows) and gene name (columns)
# Nikhil's subtypes
subtypes <- readRDS("data/milind2019/rosmap_patient_subtypes.RDS") %>%
arrange(Patient)
rownames(subtypes) <- subtypes$Patient
# Gene expression data; imputed version of syn8456719
rna <- read_tsv("data/rna_ROSMAP_DLPFC_netResidualExpression_imputed.tsv") %>%
column_to_rownames("ensembl_gene_id") %>%
dplyr::select(as.character(subtypes$Patient)) %>%
t()
# Clinical (covariates) data syn3191087
meta.cl <- read_csv("data/metadata/ROSMAP_clinical_2019-05_v3.csv") %>%
mutate(projid=as.character(projid))
# syn3382527 to match patients id's
key <- read_csv("data/metadata/ROSMAP_IDkey.csv", col_types=cols(.default = "c")) %>%
dplyr::select(projid, rnaseq_id, wgs_id) %>%
filter(rnaseq_id %in% rownames(rna)) %>%
filter(!duplicated(projid))
meta <- merge(key, meta.cl, by = "projid", all=FALSE) %>%
merge(subtypes, by.x="rnaseq_id", by.y="Patient") %>% # adding Nikil's subtyppes info
mutate(sex=if_else(msex==0, "F", "M")) %>%
dplyr::select(-(Study:spanish), -(age_at_visit_max:pmi), -individualID) %>%
arrange(rnaseq_id) # same order as data
# combined diagnosis based on syn8456629
# Remove dementia that isn't MCI or AD (cogdx=6) so that "Other" includes only MCI
meta <- meta %>%
filter(cogdx!=6) %>%
mutate(diagnosis = if_else(cogdx == 4 & braaksc >= 4 & ceradsc <= 2, "AD",
if_else(cogdx == 1 & braaksc <= 3 & ceradsc >= 3, "Control", "MCI") ))
# Same as diagnosis from Nikil's work (syn11024258)
# meta %>% dplyr::filter(diagnosis=="AD") %>% dplyr::select(Subtype) %>% unique()
# meta %>% dplyr::filter(diagnosis=="Control") %>% dplyr::select(Subtype) %>% unique()
data <- rna[meta$rnaseq_id, order(colnames(rna))]
# Translating ensembl ID to HGNC ID
ensmart <- biomaRt::useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl", verbose=TRUE)## Attempting web service request:
## http://www.ensembl.org:80/biomart/martservice?type=version&requestid=biomaRt&mart=ENSEMBL_MART_ENSEMBL
## V1
## 1 0.7
## BioMartServer running BioMart version: 0.7
## Mart virtual schema: default
## Mart host: http://www.ensembl.org:80/biomart/martservice
## Attempting web service request:
## http://www.ensembl.org:80/biomart/martservice?type=attributes&dataset=hsapiens_gene_ensembl&requestid=biomaRt&mart=ENSEMBL_MART_ENSEMBL&virtualSchema=default
## Attempting web service request:
## http://www.ensembl.org:80/biomart/martservice?type=filters&dataset=hsapiens_gene_ensembl&requestid=biomaRt&mart=ENSEMBL_MART_ENSEMBL&virtualSchema=defaultG_list<- biomaRt:::getBM(attributes=c('ensembl_gene_id', 'hgnc_symbol'), filters="ensembl_gene_id", values=colnames(data), mart=ensmart, uniqueRows = TRUE) %>%
select(Ensembl.ID=ensembl_gene_id, HGNC.ID=hgnc_symbol) %>%
mutate(HGNC.ID=if_else(HGNC.ID=="", NA_character_, HGNC.ID)) %>%
arrange(Ensembl.ID)
# Labels for plotting below
# samples assigned to a subtype are labeled as their subtype, otherwise the label is diagnosis
labels <- as.character(meta$diagnosis)
names(labels) <- meta$rnaseq_id
labels[meta$Subtype=="A"] <- "A"
labels[meta$Subtype=="B"] <- "B"PCA for all samples
There is complete overlap between AD, Control, and MCI.
It takes ~400 PC’s to cover 80% of the variance.
The first ~20 PC’s are above random based on the broken stick model, covering together 17% of the variance
pc <- prcomp(data)
if (min(pc$x[,1]) < 60) pc$x <- -pc$x
p1 <- vegan:::screeplot.prcomp(pc, bstick=TRUE, npcs = 35, main=NULL) # First ~20 PC's are above random
rel.ev <- pc$sdev/sum(pc$sdev) # proportion of variance explained
e8 <- which(cumsum(rel.ev)>0.8)[1] # it takes 395 PC's to cover 80% of the variance
plot(cumsum(rel.ev), ylab="Cumulative proportion of explained variance")
segments(0,0.8,e8,0.8, col="red", lty=4)
segments(e8,0.8,e8,0, col="red", lty=4)
autoplot(pc, data = meta, colour = 'diagnosis')
varpc <- round(100 * pc$sdev^2 / sum(pc$sdev^2), 2)Comparing polytopes
Fitting polytopes with a different number of vertices (from k=2 to k=6) using Principal Convex Hull Algorithm.
Comparing different numbre of PC’s (p=20 covering 17% and p=400 covering 80% of the variance), with and without MCI and Control samples.
200 bootstraps were generated for each space, and archetype position recalculated, indicated in the scatterplots as empty (green) dots.
For subsequent analyses, archetype position is calculated as the mean of the bootstrapped scores around each archetype.
arcfit_k3 <- list() # populated below with every combination of p and cases where k=3AD only
20 PC’s (17%)
p=20
cases <- meta$rnaseq_id[meta$diagnosis=="AD"]
pcs4arch <- t(pc$x[cases,1:p])Variance explained by different polytopes:
arc_ks = k_fit_pch(pcs4arch, ks = 2:6, check_installed = T,
bootstrap = T, bootstrap_N = 200, maxiter = 1000,
bootstrap_type = "s", seed = 2543,
volume_ratio = "none", # set to "none" if too slow
delta=0, conv_crit = 1e-04, order_type = "align",
sample_prop = 0.75)
p1 <- plot_arc_var(arc_ks, type = "varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p2 <- plot_arc_var(arc_ks, type = "res_varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p3 <- plot_arc_var(arc_ks, type = "total_var", point_size = 2, line_size = 1.5) +
theme_bw() +
ylab("Mean variance in position of vertices") # look for the highest k that gives reasonably low variance
p1
p2
p3
rm(p1, p2, p3, arc_ks)Fitting a line (k=2)
k=2
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
pl3 = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
)
rm("arcfit")Fitting a triangle (k=3)
k=3
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
p_pca = print(plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)"))
)
pl3 <- plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
)
arcfit_k3[["AD"]][["p20"]] <- arcfit
rm("arcfit")Fitting a rectangle (k=4)
k=4
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit", "pcs4arch")400 PC’s (80%)
p=400
cases <- meta$rnaseq_id[meta$diagnosis == "AD"]
pcs4arch <- t(pc$x[cases,1:p])Variance explained by different polytopes:
arc_ks = k_fit_pch(pcs4arch, ks = 2:6, check_installed = T,
bootstrap = T, bootstrap_N = 200, maxiter = 1000,
bootstrap_type = "s", seed = 2543,
volume_ratio = "none", # set to "none" if too slow
delta=0, conv_crit = 1e-04, order_type = "align",
sample_prop = 0.75)
p1 <- plot_arc_var(arc_ks, type = "varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p2 <- plot_arc_var(arc_ks, type = "res_varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p3 <- plot_arc_var(arc_ks, type = "total_var", point_size = 2, line_size = 1.5) +
theme_bw() +
ylab("Mean variance in position of vertices") # look for the highest k that gives reasonably low variance
p1
p2
p3
rm(list=c("p1", "p2", "p3", "arc_ks"))Fitting a line (k=2)
k=2
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit")Fitting a triangle (k=3)
k=3
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)"))
p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
arcfit_k3[["AD"]][["p400"]] <- arcfit
rm("arcfit")Fitting a rectangle (k=4)
k=4
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit", "pcs4arch")AD + MCI
20 PC’s (17%)
p=20
cases <- meta$rnaseq_id[meta$diagnosis %in% c("AD", "MCI")]
pcs4arch <- t(pc$x[cases,1:p])Variance explained by different polytopes:
arc_ks = k_fit_pch(pcs4arch, ks = 2:6, check_installed = T,
bootstrap = T, bootstrap_N = 200, maxiter = 1000,
bootstrap_type = "s", seed = 2543,
volume_ratio = "none", # set to "none" if too slow
delta=0, conv_crit = 1e-04, order_type = "align",
sample_prop = 0.75)
p1 <- plot_arc_var(arc_ks, type = "varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p2 <- plot_arc_var(arc_ks, type = "res_varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p3 <- plot_arc_var(arc_ks, type = "total_var", point_size = 2, line_size = 1.5) +
theme_bw() +
ylab("Mean variance in position of vertices") # look for the highest k that gives reasonably low variance
p1
p2
p3
rm(list=c("p1", "p2", "p3", "arc_ks"))Fitting a line (k=2)
k=2
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit")Fitting a triangle (k=3)
k=3
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)"))
p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
arcfit_k3[["ADMCI"]][["p20"]] <- arcfit
rm("arcfit")Fitting a rectangle (k=4)
k=4
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit", "pcs4arch")400 PC’s (80%)
p=400
cases <- meta$rnaseq_id[meta$diagnosis %in% c("AD", "MCI")]
pcs4arch <- t(pc$x[cases,1:p])Variance explained by different polytopes:
arc_ks = k_fit_pch(pcs4arch, ks = 2:6, check_installed = T,
bootstrap = T, bootstrap_N = 200, maxiter = 1000,
bootstrap_type = "s", seed = 2543,
volume_ratio = "none", # set to "none" if too slow
delta=0, conv_crit = 1e-04, order_type = "align",
sample_prop = 0.75)
p1 <- plot_arc_var(arc_ks, type = "varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p2 <- plot_arc_var(arc_ks, type = "res_varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p3 <- plot_arc_var(arc_ks, type = "total_var", point_size = 2, line_size = 1.5) +
theme_bw() +
ylab("Mean variance in position of vertices") # look for the highest k that gives reasonably low variance
p1
p2
p3
rm(list=c("p1", "p2", "p3", "arc_ks"))Fitting a line (k=2)
k=2
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit")Fitting a triangle (k=3)
k=3
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)"))
p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
arcfit_k3[["ADMCI"]][["p400"]] <- arcfit
rm("arcfit")Fitting a rectangle (k=4)
k=4
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit", "pcs4arch")All (AD + MCI + Control)
20 PC’s (17%)
p=20
cases <- meta$rnaseq_id[meta$diagnosis %in% c("AD", "MCI", "Control")]
pcs4arch <- t(pc$x[cases,1:p])Variance explained by different polytopes:
arc_ks = k_fit_pch(pcs4arch, ks = 2:6, check_installed = T,
bootstrap = T, bootstrap_N = 200, maxiter = 1000,
bootstrap_type = "s", seed = 2543,
volume_ratio = "none", # set to "none" if too slow
delta=0, conv_crit = 1e-04, order_type = "align",
sample_prop = 0.75)
p1 <- plot_arc_var(arc_ks, type = "varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p2 <- plot_arc_var(arc_ks, type = "res_varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p3 <- plot_arc_var(arc_ks, type = "total_var", point_size = 2, line_size = 1.5) +
theme_bw() +
ylab("Mean variance in position of vertices") # look for the highest k that gives reasonably low variance
p1
p2
p3
rm(list=c("p1", "p2", "p3", "arc_ks"))Fitting a line (k=2)
k=2
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit")Fitting a triangle (k=3)
k=3
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)"))
p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
arcfit_k3[["All"]][["p20"]] <- arcfit
rm("arcfit")Fitting a rectangle (k=4)
k=4
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit", "pcs4arch")400 PC’s (80%)
p=400
cases <- meta$rnaseq_id[meta$diagnosis %in% c("AD", "MCI", "Control")]
pcs4arch <- t(pc$x[cases,1:p])Variance explained by different polytopes:
arc_ks = k_fit_pch(pcs4arch, ks = 2:6, check_installed = T,
bootstrap = T, bootstrap_N = 200, maxiter = 1000,
bootstrap_type = "s", seed = 2543,
volume_ratio = "none", # set to "none" if too slow
delta=0, conv_crit = 1e-04, order_type = "align",
sample_prop = 0.75)
p1 <- plot_arc_var(arc_ks, type = "varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p2 <- plot_arc_var(arc_ks, type = "res_varexpl", point_size = 2, line_size = 1.5) + theme_bw()
p3 <- plot_arc_var(arc_ks, type = "total_var", point_size = 2, line_size = 1.5) +
theme_bw() +
ylab("Mean variance in position of vertices") # look for the highest k that gives reasonably low variance
p1
p2
p3
rm(list=c("p1", "p2", "p3", "arc_ks"))Fitting a line (k=2)
k=2
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit")Fitting a triangle (k=3)
k=3
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)"))
p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
arcfit_k3[["All"]][["p400"]] <- arcfit
rm("arcfit")Fitting a rectangle (k=4)
k=4
arcfit <- fit_pch_bootstrap(pcs4arch, n = 200, sample_prop = 0.75, seed = 2543,
noc = k, delta = 0, conv_crit = 1e-04, type = "s")
# empty points are bootstrapped data showing variance around each archetype
plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:3, line_size = 1.5,
text_size = 24, data_size = 3,
data_lab = labels[cases]
#colors= palette(rainbow(6))
) p_pca = plot_arc(arc_data = arcfit, data = pcs4arch,
which_dimensions = 1:2, line_size = 1.5,
data_lab = labels[cases],
text_size = 6, data_size = 3
)
rm("arcfit", "pcs4arch")Distances and classification k=3
Calculating Euclidean distance of each sample (including AD, Control, and MCI) from each of the three archetypes.
Each sample is classified to the nearest archetype.
k=3
archDistClass <- list()
for (li in 1:length(arcfit_k3)) {
for (lj in 1:length(arcfit_k3[[li]])) {
arcfit <- arcfit_k3[[li]][[lj]]
# position of archetypes as the mean of the bootstrapped samples:
Xpca <- t(average_pch_fits(arcfit)$XC)
rownames(Xpca) <- c("Archetype_1", "Archetype_2", "Archetype_3")
allscores <- pc$x[,1:ncol(Xpca)]
# euclidean distances for each sample from each archetype
archetypes <- as.matrix(dist(rbind(Xpca, allscores)), method="euclidean")[-c(1:3),1:3] %>%
as.data.frame()
# classify each sample to closest archetype
archetypes$Archetype <- apply(archetypes, 1, function(x){names(x)[x==min(x)]})
archetypes$rnaseq_id <- rownames(archetypes)
nli <- names(arcfit_k3)[[li]]
nlj <- names(arcfit_k3[[li]])[[lj]]
archDistClass[[nli]][[nlj]] <- archetypes
if (li==1 & lj==1) {arch.pca <- rbind(Xpca, allscores)}
rm(archetypes)
}
}df <- as.data.frame(arch.pca[-c(1:3),])
df$Archetype <- archDistClass[["AD"]][["p20"]]$Archetype
df$Diagnosis <- meta$diagnosis
archpos <- as.data.frame(arch.pca[1:3,]) %>% rownames_to_column(var="Archetype")
ggplot(df) +
geom_point(aes(x=PC1, y=PC2, color=Archetype, shape=Diagnosis)) +
geom_point(data=archpos, aes(x=PC1, y=PC2), color="black", size=10, shape="*") +
geom_label(data=archpos, aes(x=PC1, y=PC2, label=c("1","2","3")), nudge_y = 3, size=3) +
xlab(paste0("PC1 (", varpc[1], "%)")) +
ylab(paste0("PC2 (", varpc[2], "%)")) +
ggtitle("Archetype classfication and position; k=3 p=20 AD only")
Comparing distances
Comparing distances from each archetype for the different combinations of p and cases, when k=3
Archetype 1
arch <- "Archetype_1"
Adist <- c()
for (anm in names(archDistClass)) {
for (pnm in names(archDistClass[[1]])) {
grnm <- paste(anm, pnm, sep="_")
Adist[[grnm]] <- archDistClass[[anm]][[pnm]][[arch]]
}
}
ggpairs(as.data.frame(Adist), title="Comparing distances - Archetype 1") 
Archetype 2
arch <- "Archetype_2"
Adist <- c()
for (anm in names(archDistClass)) {
for (pnm in names(archDistClass[[1]])) {
grnm <- paste(anm, pnm, sep="_")
Adist[[grnm]] <- archDistClass[[anm]][[pnm]][[arch]]
}
}
ggpairs(as.data.frame(Adist), title="Comparing distances - Archetype 2") 
Archetype 3
arch <- "Archetype_3"
Adist <- c()
for (anm in names(archDistClass)) {
for (pnm in names(archDistClass[[1]])) {
grnm <- paste(anm, pnm, sep="_")
Adist[[grnm]] <- archDistClass[[anm]][[pnm]][[arch]]
}
}
ggpairs(as.data.frame(Adist), title="Comparing distances - Archetype 3") 
Comparing classification
Aclass <-list()
for (anm in names(archDistClass)) {
for (pnm in names(archDistClass[[1]])) {
grnm <- paste(anm, pnm, sep="_")
Aclass[[grnm]] <- archDistClass[[anm]][[pnm]][["Archetype"]]
}
}
ggpairs(as.data.frame(Aclass), title="Comparing classification", lower=list(discrete="ratio")) 
Top genes k=3 p=20 AD only
Get the genes significantly associated with each archetype for k=3, p=20, and AD cases only.
Significance is determined by fitting a decreasing function with ditance from a given archetype to the expression level of each gene. Therefore, the set of top genes associated with each archetypes are genes that are significantly upregulated in that archetype relative to the other two archetypes.
arcfit <- arcfit_k3[["AD"]][["p20"]]
cases <- meta$rnaseq_id[meta$diagnosis %in% c("AD", "MCI", "Control")]
pcs4arch <- t(pc$x[cases,1:20])
data_attr <- average_pch_fits(arcfit) %>%
merge_arch_dist(data = pcs4arch, feature_data = t(data),
dist_metric = "euclidean", rank = F)
enriched_genes = find_decreasing_wilcox(data_attr$data, data_attr$arc_col,
features = data_attr$features_col,
bin_prop = 0.1, method = "BioQC")
# get genes that are a decreasing function of distance from either archetypes
# p < 0.01
topgenes = get_top_decreasing(summary_genes = enriched_genes,
cutoff_genes = 0.01, cutoff_sets = 0.05,
cutoff_metric = "wilcoxon_p_val",
p.adjust.method = "fdr",
order_by = "mean_diff", order_decreasing = T,
min_max_diff_cutoff_g = 0.4, min_max_diff_cutoff_f = 0.03)## -- archetype_3
##
## ENSG00000147571, ENSG00000128564, ENSG00000186212, ENSG00000183090
## ENSG00000263667, ENSG00000197322, ENSG00000120875, ENSG00000086717
## ENSG00000150175, ENSG00000164600, ENSG00000141433, ENSG00000136750
##
## -- archetype_1
##
## ENSG00000159167, ENSG00000091879, ENSG00000196154, ENSG00000010310
## ENSG00000204389, ENSG00000108691, ENSG00000228058, ENSG00000152049
## ENSG00000275395, ENSG00000173530, ENSG00000184557, ENSG00000106211
##
## -- archetype_2
##
## ENSG00000244734, ENSG00000134817, ENSG00000047457, ENSG00000152661
## ENSG00000141469, ENSG00000110436, ENSG00000171885, ENSG00000244731
## ENSG00000265972, ENSG00000224389, ENSG00000164089, ENSG00000162407topgenes <- topgenes[["enriched_genes"]] %>%
rename(Archetype=arch_name, Ensembl.ID=genes) %>%
mutate(Cohort="rosmap", Archetype=str_to_sentence(Archetype)) %>%
left_join(G_list, by="Ensembl.ID")
table(topgenes$Archetype)##
## Archetype_1 Archetype_2 Archetype_3
## 519 463 937# overlap of topgenes between archetypes
knitr::kable(filter(topgenes, Archetype!="Archetype_3") %>%
filter(duplicated(Ensembl.ID)) %>%
dplyr::select(Ensembl.ID, HGNC.ID),
caption = "Overlap between archetype_1 and archetype_2:")| Ensembl.ID | HGNC.ID |
|---|---|
| ENSG00000126803 | HSPA2 |
| ENSG00000134569 | LRP4 |
| ENSG00000168209 | DDIT4 |
| ENSG00000164188 | RANBP3L |
| ENSG00000117318 | ID3 |
| ENSG00000174607 | UGT8 |
| ENSG00000181826 | RELL1 |
| ENSG00000198959 | TGM2 |
| ENSG00000102760 | RGCC |
| ENSG00000158825 | CDA |
| ENSG00000250722 | SELENOP |
| ENSG00000177807 | KCNJ10 |
| ENSG00000117594 | HSD11B1 |
| ENSG00000169715 | MT1E |
| ENSG00000102362 | SYTL4 |
| ENSG00000122862 | SRGN |
| ENSG00000104267 | CA2 |
knitr::kable(filter(topgenes, Archetype!="Archetype_2") %>%
filter(duplicated(Ensembl.ID)) %>%
dplyr::select(Ensembl.ID, HGNC.ID),
caption = "Overlap between archetype_1 and archetype_3:")| Ensembl.ID | HGNC.ID |
|---|
knitr::kable(filter(topgenes, Archetype!="Archetype_1") %>%
filter(duplicated(Ensembl.ID)) %>%
dplyr::select(Ensembl.ID, HGNC.ID),
caption = "Overlap between archetype_3 and archetype_2:")| Ensembl.ID | HGNC.ID |
|---|---|
| ENSG00000197921 | HES5 |
logFC k=3 p=20 AD only
A sample of “controls” is generated based on the centroid of the triangle formed by the three archetypes in each bootstrapped sapce. This results in a total of 199 “control” replicates in addition to the 199 replicates for each archetype (4 “populations” in total).
Everything is then rotated back to the original gene expression space, and limma is applied to each of the archetyppe population relative to the control population.
The rationale behind it is that the top genes are defined as those that are upregulated for a given archetype relative to the other two, implying that the relative baseline is around the mean of the three archetypes.
arcfit <- arcfit_k3[["AD"]][["p20"]]
k = ncol(arcfit$pch_fits$XC[[1]]) # number of archetypes
p = nrow(arcfit$pch_fits$XC[[1]]) # number of PC's
nb = length(arcfit$pch_fits$XC) # number of bootstrapped replicates
archnames = c("archetype_1", "archetype_2", "archetype_3")
# calculate PC scores for the centroid ("control" samples) in each bootstrapped space
Xcontrol <- arcfit$pch_fits$XC %>%
sapply(rowMeans) %>%
t()
# Reorganize bootstrapped scores into a list of 4 nb x p matrices
# one matrix for each archetype and one for control
XL <- arcfit$pch_fits$XC %>%
unlist() %>%
array(dim=c(p,k,nb)) %>%
aperm(perm=c(3,1,2)) %>%
abind::abind(Xcontrol) %>%
plyr::alply(3,.dims = TRUE)
names(XL) <- c(archnames, "Control")
# rotate everything back to original space
# and translate to original center
# (reversing the pca procedure)
rotM <- pc$rotation # rotation matrix (eigenvectors)
mu <- pc$center # means of gene expression data by which pc's were centered
boot.ge <- list()
for (pop in names(XL)) {
Xbp <- XL[[pop]]
boot.ge[[pop]] <- Xbp %*% t(rotM[,1:p]) %>% # rotating
scale(center = -mu, scale = FALSE) # translating
rm(Xbp)
}
# calculate LogFC using limma
lm.res <- list()
for (a in archnames) {
type <- as.factor(rep(c(a,"control"), each=nb)) # disease status
de.df <- as.data.frame(t(rbind(boot.ge[[a]], boot.ge[["Control"]])))
fit <- limma::lmFit(de.df, design=model.matrix(~type))
fit <- limma::eBayes(fit)
tt <- limma::topTable(fit, number=Inf, coef=2) %>%
rownames_to_column("Ensembl.ID") %>%
left_join(G_list)
lm.res[[a]] <- tt
rm(tt,fit, de.df, type)
}
logFC <- data.frame(Ensembl.ID = lm.res[[1]]$Ensembl.ID,
Archetype_1 = lm.res[[1]]$logFC,
Archetype_2 = lm.res[[2]]$logFC,
Archetype_3 = lm.res[[3]]$logFC)
# visualize with heatmap
ComplexHeatmap::Heatmap(as.matrix(logFC[,-1]),
show_row_dend = FALSE,
show_row_names = FALSE,
)
Output
Merge logFC with topgenes keeping logFC for only topgenes
topgenes <- left_join(topgenes,
logFC %>% pivot_longer(-Ensembl.ID, names_to = "Archetype", values_to = "logFC"),
by=c("Ensembl.ID", "Archetype"),
all.x=TRUE, all.y=FALSE)Merge gene expression profile for each archetype with the original gene expression dataset
arch.ge <- t(sapply(boot.ge, colMeans)[,1:3]) %>%
rbind(data)Merge archetype distances and classification with metadata
archetypes_meta <- left_join(meta, archDistClass[["AD"]][["p20"]], by="rnaseq_id", all=TRUE) %>%
rename(Cognitive.Diagnosis=cogdx, Sex=sex, Diagnosis=diagnosis, Braak.Score=braaksc, CERAD.Score=ceradsc, APOE4=apoe_genotype)Save objects
outpath <- "analyses/rnaseq/1_fit_archetypes/"
saveRDS(arcfit_k3, file=paste0(outpath, "arcfit_k3.RDS"))
saveRDS(archetypes_meta, file=paste0(outpath, "archetypes_meta_k3p20_AD.RDS"))
saveRDS(list(arch.pca=arch.pca, all.pca=pc), file=paste0(outpath, "archetypes_pca_k3p20_AD.RDS"))
saveRDS(arch.ge, file=paste0(outpath, "archetypes_ge_k3p20_AD.RDS"))
saveRDS(topgenes, file=paste0(outpath, "topgenes_logFC_k3p20_AD.RDS"))
saveRDS(lm.res, file=paste0(outpath, "logFC_allres_k3p20_AD.RDS"))Session info
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 3.6.2 (2019-12-12)
## os macOS Mojave 10.14.6
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2020-09-09
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib
## AnnotationDbi 1.48.0 2019-10-29 [1]
## AnnotationHub 2.18.0 2019-10-29 [1]
## assertthat 0.2.1 2019-03-21 [1]
## backports 1.1.6 2020-04-05 [1]
## Biobase 2.46.0 2019-10-29 [1]
## BiocFileCache 1.10.2 2019-11-08 [1]
## BiocGenerics 0.32.0 2019-10-29 [1]
## BiocManager 1.30.10 2019-11-16 [1]
## BiocVersion 3.10.1 2019-06-06 [1]
## BioQC 1.14.0 2019-10-29 [1]
## bit 1.1-15.2 2020-02-10 [1]
## bit64 0.9-7 2017-05-08 [1]
## blob 1.2.1 2020-01-20 [1]
## broom 0.5.5 2020-02-29 [1]
## callr 3.4.3 2020-03-28 [1]
## cellranger 1.1.0 2016-07-27 [1]
## cli 2.0.2 2020-02-28 [1]
## codetools 0.2-16 2018-12-24 [1]
## colorspace 1.4-1 2019-03-18 [1]
## cowplot * 1.0.0 2019-07-11 [1]
## crayon 1.3.4 2017-09-16 [1]
## curl 4.3 2019-12-02 [1]
## data.table * 1.12.8 2019-12-09 [1]
## DBI 1.1.0 2019-12-15 [1]
## dbplyr 1.4.2 2019-06-17 [1]
## desc 1.2.0 2018-05-01 [1]
## devtools 2.3.0 2020-04-10 [1]
## digest 0.6.25 2020-02-23 [1]
## dplyr * 0.8.5 2020-03-07 [1]
## edgeR 3.28.1 2020-02-26 [1]
## ellipsis 0.3.0 2019-09-20 [1]
## evaluate 0.14 2019-05-28 [1]
## fansi 0.4.1 2020-01-08 [1]
## fastmap 1.0.1 2019-10-08 [1]
## forcats * 0.5.0 2020-03-01 [1]
## fs 1.4.1 2020-04-04 [1]
## generics 0.0.2 2018-11-29 [1]
## GGally * 1.5.0 2020-03-25 [1]
## ggfortify * 0.4.9 2020-03-11 [1]
## ggplot2 * 3.3.0 2020-03-05 [1]
## glue 1.4.0 2020-04-03 [1]
## gridExtra 2.3 2017-09-09 [1]
## gtable 0.3.0 2019-03-25 [1]
## haven 2.2.0 2019-11-08 [1]
## highr 0.8 2019-03-20 [1]
## hms 0.5.3 2020-01-08 [1]
## htmltools 0.4.0 2019-10-04 [1]
## htmlwidgets 1.5.1 2019-10-08 [1]
## httpuv 1.5.2 2019-09-11 [1]
## httr 1.4.1 2019-08-05 [1]
## interactiveDisplayBase 1.24.0 2019-10-29 [1]
## IRanges 2.20.2 2020-01-13 [1]
## jsonlite 1.6.1 2020-02-02 [1]
## knitr 1.28 2020-02-06 [1]
## later 1.0.0 2019-10-04 [1]
## lattice 0.20-41 2020-04-02 [1]
## lazyeval 0.2.2 2019-03-15 [1]
## lifecycle 0.2.0 2020-03-06 [1]
## limma 3.42.2 2020-02-03 [1]
## locfit 1.5-9.4 2020-03-25 [1]
## lpSolve * 5.6.15 2020-01-24 [1]
## lubridate 1.7.8 2020-04-06 [1]
## magrittr 1.5 2014-11-22 [1]
## Matrix * 1.2-18 2019-11-27 [1]
## memoise 1.1.0 2017-04-21 [1]
## mime 0.9 2020-02-04 [1]
## modelr 0.1.6 2020-02-22 [1]
## munsell 0.5.0 2018-06-12 [1]
## nlme 3.1-145 2020-03-04 [1]
## ParetoTI * 0.1.13 2020-04-11 [1]
## pillar 1.4.3 2019-12-20 [1]
## pkgbuild 1.0.6 2019-10-09 [1]
## pkgconfig 2.0.3 2019-09-22 [1]
## pkgload 1.0.2 2018-10-29 [1]
## plotly 4.9.2.1 2020-04-04 [1]
## plyr 1.8.6 2020-03-03 [1]
## prettyunits 1.1.1 2020-01-24 [1]
## processx 3.4.2 2020-02-09 [1]
## promises 1.1.0 2019-10-04 [1]
## ps 1.3.2 2020-02-13 [1]
## purrr * 0.3.3 2019-10-18 [1]
## R6 2.4.1 2019-11-12 [1]
## rappdirs 0.3.1 2016-03-28 [1]
## RColorBrewer 1.1-2 2014-12-07 [1]
## Rcpp 1.0.4.6 2020-04-09 [1]
## readr * 1.3.1 2018-12-21 [1]
## readxl 1.3.1 2019-03-13 [1]
## remotes 2.1.1 2020-02-15 [1]
## reprex 0.3.0 2019-05-16 [1]
## reshape 0.8.8 2018-10-23 [1]
## reticulate * 1.15 2020-04-02 [1]
## rlang 0.4.5 2020-03-01 [1]
## rmarkdown 2.1 2020-01-20 [1]
## rprojroot 1.3-2 2018-01-03 [1]
## RSQLite 2.2.0 2020-01-07 [1]
## rstudioapi 0.11 2020-02-07 [1]
## rvest 0.3.5 2019-11-08 [1]
## S4Vectors 0.24.4 2020-04-09 [1]
## scales 1.1.0 2019-11-18 [1]
## sessioninfo 1.1.1 2018-11-05 [1]
## shiny 1.4.0.2 2020-03-13 [1]
## stringi 1.4.6 2020-02-17 [1]
## stringr * 1.4.0 2019-02-10 [1]
## testthat 2.3.2 2020-03-02 [1]
## tibble * 3.0.0 2020-03-30 [1]
## tidyr * 1.0.2 2020-01-24 [1]
## tidyselect 1.0.0 2020-01-27 [1]
## tidyverse * 1.3.0 2019-11-21 [1]
## usethis 1.6.0 2020-04-09 [1]
## vctrs 0.2.4 2020-03-10 [1]
## viridisLite 0.3.0 2018-02-01 [1]
## withr 2.1.2 2018-03-15 [1]
## xfun 0.12 2020-01-13 [1]
## xml2 1.3.1 2020-04-09 [1]
## xtable 1.8-4 2019-04-21 [1]
## yaml 2.2.1 2020-02-01 [1]
## source
## Bioconductor
## Bioconductor
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## Bioconductor
## Bioconductor
## Bioconductor
## CRAN (R 3.6.0)
## Bioconductor
## Bioconductor
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.1)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.1)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## Bioconductor
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## Bioconductor
## Bioconductor
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## Bioconductor
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## Github (vitkl/ParetoTI@5109906)
## CRAN (R 3.6.1)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.1)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## Bioconductor
## CRAN (R 3.6.1)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.1)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.2)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
##
## [1] /Users/habera/Dropbox/ampad_archetypes/renv/library/R-3.6/x86_64-apple-darwin15.6.0
## [2] /private/var/folders/ls/45q5krvs5xd4m27ks7cmhy80g4hpwp/T/Rtmpa30vBk/renv-system-library
## [3] /Library/Frameworks/R.framework/Versions/3.6/Resources/library